%% Matrix multiplication =====================================================
clc
A = randn(100, 100);
B = randn(100, 100);
% using built-in function
tic
C = A*B;
toc
% is more efficient than
tic
C = zeros(100, 100);
for i = 1:100
    for j = 1:100
        for k = 1:100
            C(i,j) = C(i,j) + A(i,k) * B(k,j);
        end
    end
end
toc

% Why? Because when you use the proper symbol for the operation
% you mean to compute, Matlab knows what you're doing and it
% can call highly optimized routines in C or Fortran that
% have been engineered for decades to takes advantages
% of your CPU architecture.
% Part of this is called "SIMD": single instruction, multiple data.
% Think about it: when multiplying two matrices, you are doing the
% same operations on many different numbers. CPUs can do many of
% those in parallel, /if/ they know that this is what you want.
%
% Read more here:
% https://ch.mathworks.com/help/matlab/matlab_prog/vectorization.html


%% Multiplying each column of a matrix by a different scalar ====================
clc
A = randn(100, 100);
x = randn(100, 1);
% using matrix product 
tic
B = A*diag(x);
toc
% is more efficient than
tic
B = zeros(100, 100);
for i = 1:100
    B(:,i) = A(:,i)*x(i);
end
toc
% BUT: see also "broadcasting.m" for even faster code.


%% Unary function applied to a vector =================================
clc
x = randn(100, 1);
% using vectorized function
tic
y = sin(x); % applying the same function to many numbers
toc
% is more efficient than
tic
y = zeros(100, 1);
for i = 1:100
    y(i) = sin(x(i));
end
toc

%% Computing aggregate statistics ============================================
clc
x = randn(100, 1);
% using vectorized functions
tic
n = length(x);
mu = sum(x)/n;
sigma = sqrt(sum((x-mu).^2)/(n-1));
toc
% is more efficient than
tic
n = length(x);
mu = 0;
sigma = 0;
for i = 1:n
    mu = mu + x(i);
end
mu = mu/n;
for i = 1:n
    sigma = sigma + (x(i)-mu)^2;
end
sigma = sqrt(sigma/(n-1));
toc

% Note: here, you can also use mean(x) and std(x).

%% Compute cumulative product with a condition ==========================
clc
n = 10;
x = rand(1, n);

% Naive approach (using a loop)
y_naive = ones(1, n);

tic
for i = 2:n
    if x(i) > 0.5
        y_naive(i) = y_naive(i-1) * x(i);
    else
        y_naive(i) = y_naive(i-1);
    end
end
toc

% Vectorized approach with well-chosen indexing
tic
mask = (x > 0.5); % take a look inside!
y_vectorized = cumprod([1, x(mask)]);
y_vectorized = y_vectorized(cumsum(mask) + 1);
toc